home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / igamma.pro < prev    next >
Encoding:
Text File  |  1997-07-08  |  3.6 KB  |  117 lines
  1. ;$Id: igamma.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       IGAMMA
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the incomplete gamma function, Px(a).
  11. ;
  12. ; CATEGORY:
  13. ;       Special Functions.
  14. ;
  15. ; CALLING SEQUENCE:
  16. ;       Result = Igamma(a, x)
  17. ;
  18. ; INPUTS:
  19. ;       A:    A positive scalar of type integer, float or double that
  20. ;             specifies the parametric exponent of the integrand.
  21. ;
  22. ;       X:    A positive scalar of type integer, float or double that 
  23. ;             specifies the upper limit of integration.
  24. ;
  25. ; KEYWORD PARAMETERS:
  26. ;       METHOD:  Use this keyword to specify a named variable which returns
  27. ;                the method used to compute the incomplete gamma function.
  28. ;                A value of 0 indicates that a power series representation
  29. ;                was used. A value of 1 indicates that a continued fractions
  30. ;                method was used.
  31. ;
  32. ; EXAMPLE:
  33. ;       Compute the incomplete gamma function for the corresponding elements
  34. ;       of A and X.
  35. ;       Define the parametric exponents.
  36. ;         A = [0.10, 0.50, 1.00, 1.10, 6.00, 26.00]
  37. ;       Define the the upper limits of integration.
  38. ;         X = [0.0316228, 0.0707107, 5.00000, 1.04881, 2.44949, 25.4951]
  39. ;       Allocate an array to store the results.
  40. ;         result = fltarr(n_elements(A))
  41. ;       Compute the incomplete gamma functions.
  42. ;         for k = 0, n_elements(A)-1 do $
  43. ;           result(k) = Igamma(A(k), X(k))
  44. ;       The result should be:
  45. ;         [0.742026, 0.293128, 0.993262, 0.607646, 0.0387318, 0.486387]
  46. ;
  47. ; PROCEDURE:
  48. ;       IGAMMA computes the incomplete gamma function, Px(a), using either
  49. ;       a power series representation or a continued fractions method. If X
  50. ;       is less than or equal to A+1, a power series representation is used.
  51. ;       If X is greater than A+1, a continued fractions method is used. 
  52. ; 
  53. ; REFERENCE:
  54. ;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
  55. ;       Cambridge University Press
  56. ;       ISBN 0-521-43108-5
  57. ;
  58. ; MODIFICATION HISTORY:
  59. ;       Written by:  GGS, RSI, September 1994
  60. ;                    IGAMMA is based on the routines: gser.c, gcf.c, and  
  61. ;                    gammln.c described in section 6.2 of Numerical Recipes,
  62. ;                    The Art of Scientific Computing (Second Edition), and is
  63. ;                    used by permission.
  64. ;       Modified:    GGS, RSI, January 1996
  65. ;                    Corrected the case of IGAMMA(a, 0) for a > 0.
  66. ;-
  67.  
  68. function igamma, a, x, itmax = itmax, method = method
  69.  
  70.   on_error, 2
  71.  
  72.   if a le 0 or x lt 0 then $
  73.     message, 'a and x must be positive scalars.'
  74.  
  75.   if x eq 0.0 then return, 0.0
  76.  
  77.   eps = 3.0e-7
  78.   fpmin = 1.0e-30
  79.   if keyword_set(itmax) eq 0 then itmax = 100
  80.  
  81.   if x le (a + 1) then begin ;Series Representation.
  82.     method = 0
  83.     ap = a
  84.     sum = 1.0 / a
  85.     del = sum
  86.     for k = 1, itmax do begin
  87.       ap = ap + 1.0
  88.       del = del * x / ap
  89.       sum = sum + del
  90.       if abs(del) lt abs(sum)*eps then return, $
  91.         sum * exp(-x + a * alog(x) - lngamma(a))
  92.     endfor 
  93.   endif else begin ;Continued Fractions.
  94.     method = 1
  95.     b = x + 1.0 - a
  96.     c = 1.0 / fpmin
  97.     d = 1.0 / b
  98.     h = d
  99.     for k = 1, itmax do begin
  100.       an = -k * (k - a)
  101.       b = b + 2
  102.       d = an * d + b
  103.       if abs(d) lt fpmin then d = fpmin
  104.       c = b + an / c
  105.       if abs(c) lt fpmin then c = fpmin
  106.       d = 1.0 / d
  107.       del = d * c
  108.       h = h * del
  109.       if abs(del - 1) lt eps then return, $
  110.         1 - (exp(-x + a * alog(x) - lngamma(a)) * h)
  111.     endfor
  112.   endelse
  113.  
  114.   message, 'Failed to converge within given parameters.'
  115.  
  116. end
  117.